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Classical color fields produced by the small- x wave functions of colliding ultrarelativistic nuclei 
ON \ have been numerically computed. We set up the framework for computing the production of 

small-mass quark-antiquark pairs in these color fields by numerically integrating the Dirac 
equation. This computation is essential for understanding the conversion of the initial gluonic 
state to chemically equilibrated quark-gluon plasma. To illustrate and overcome technical 
difficulties associated with the longitudinal dimension, we first consider numerically the case 
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of one time + one longitudinal space dimension. 
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1 Introduction 



The dynamics of an ultrarelativistic heavy ion collisions is usually described in the following 
terms: two nuclei in their T = 0, entropy=0 ground state move along the light cone, collide 
at t = and form a large-entropy extended system with deconfined quark-gluon degrees of 
freedom. This system expands, passes a QCD phase transition, converts itself to a hadronic 
phase, which finally decouples and sends hadrons to detectors. 

Important partial confirmation for this scenario comes from recent experimental results 
from the relativistic heavy ion collider RHIC. These suggest that in y/s = 200 GeV Au+Au 
collisions a nearly thermalised quark-gluon plasma is formed pQ. One of the main pieces of 
evidence comes from azimuthal asymmetries in non-central collisions [21IH1II] : a hydrodynamic 
computation [2] , with an assumed equation of state and initial conditions fitted to transverse 
spectra, shows that the initial spatial azimuthal asymmetry is converted to just the correct 
amount of momentum space azimuthal asymmetry if the equation of state is that of an ideal 
fluid. 

There is one significant deficiency in the theoretical analysis of this scenario: virtually all 
the models describing the initial state (for examples, see [HI 13 El El El) are based on the 
almost purely gluonic small- x partonic content of the nuclear wave function, while an ideal 
quark-gluon plasma would contain gluons and quarks+antiquarks in the ratio 16/(21iVy/2) 
and anyway the final hadronic state contains flavour in a fully thermalised manner At 
what stage do the small-mass u,d, and s flavour degrees of freedom appear in the system? 
Experiments do not yet shed any light on this problem. 

Models based on weakly coupled quark-gluon degrees of freedom, like parton cascade mod- 
els, fail to reproduce both kinetic and chemical equilibration J^j; the coupling is so weak 
that collision times become too large relative to the lifetime of the system. The purpose 
of this paper is to start from a strongly coupled and phenomenologically viable model, the 
classical field computation of gluon production in a collision of two nuclei 1131 114j using 
the McLerran-Venugopalan model 8 for the distribution of gluons in a single nucleus, later 
evolved and termed color glass condensate (see arid references therein) . We then discuss 
how the amount of small-mass u,d, and s quark-antiquark pairs produced by the colour fields 
of this model can be computed by numerically integrating the evolution of a negative energy 
spinor as given by the Dirac equation and projecting on a positive energy spinorfTHl I17j. In 
this paper we give numerical results only for a 1+1 dimensional toy model version of the full 
computation, to establish the viability of the method. 

The following should be emphasised from the outset: 

• This is not a computation of pair production from strong colour fields by quantum 
tunneling via the Schwinger mechanism, which has often been studied |18U19j . Instead, 
the pairs are produced via multiple interactions of quasi-real Fourier components of the 
color fields; in the dilute limit this is just the two gluon fusion mechanism g*+g*^ 
q+q , which for heavy quarks also is the dominant mechanism [20]. The same produc- 
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tion mechanism has been studied in where several approximations for the quark 
retarded propagator in an external field have been investigated. 

• Basically, there are two quantities we would like to know: how fast does the qq density 
grow in comparison with the gluon density (what are the typical production times 
in units of 1/Q S , Qs = saturation scale, see below) and how high is the qq density in 
comparison with the gluon density (what are the total energies per unit rapidity in units 
of R^Qg). Parametrically, quark pair production is suppressed by a factor a s , but we 
are not in the weak coupling limit. In the strong coupling regime, a large qq component 
could be created, as required for chemical equilibration. Kinetic equilibration of the 
longitudinal degree of freedom is still an open issue. 

• The pair production being computed in a given colour field, the feedback is not taken 
into account. The results will thus be quantitatively reliable only as long as the energy 
in qq pairs remains less than that in gluons. For the Abelian Higgs model in 1+1 
dimensions and with different initial conditions a numerical scheme for including both 
bosonic and fermionic dynamical degrees of freedom has been developed in |22| . 

• The computation of the qq production is technically much more complicated than that 
of gluons. The colour fields giving rise to gluons are assumed to be independent of the 
space-time rapidity rj = | log(x + /x~), they only depend on t,Xt- Thus strict boost 
invariance for gluons is obtained (unless rapidity dependence is introduced via that of 
the saturation scales) and the single rapidity of the problem, that of the gluon, can be 
completely removed from the equations. For quark pair production, two rapidities enter 
and a non-trivial dependence in Ay = y q — yq appears. This also implies that the quark 
wave function ip(r, rj, Xy) will depend on all the 3+1 variables. However, formulating 
the initial condition on the light cones (say, x~ = 0, x + > 0) is impossible using the 
natural variables t, 77 since fixed x + = r exp(^)/v2 > cannot be reached for r — > 
unless also rj — > 00. We thus have to use as variables r, x^ or, more symmetrically, r, z. 

• For an approach to quark pair production via special nonperturbative instanton con- 
figurations, see 

The structure of this paper is as follows. In section [21 the problem is formulated in full 
generality in 3+1 dimensions. Particular attention is given to the initial condition and the 
difficulties associated with the longitudinal dimension are pointed out. This leads us to 
truncate the full theory by neglecting all the transverse integrations to a 1+1 dimensional toy 
model, with which we can test the numerical solution of the time+longitudinal dependence. 
The free Dirac equation in 1+1 dimensions using (r, rj), (t,x^) or (r, z) as variables (we do 
not go all the way to (t, z)\) is studied and solved analytically in section |21 Its numerical 
solution is carried out in section |1] and shown to agree with the analytic one. Finally, in 
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section the 1+1 dimensional Dirac equation is solved with various forms of the external 
gluonic field. 

In the time- longitudinal space we shall use three sets of variables: t, z, ds 2 = dt 2 — dz 2 , 
the light cone coordinates x = (t ± z)/\/2 = re^/V^j ds 2 = 2dx + dx~ and proper time 
and spacetime rapidity r = \/t 2 — z 2 = \Jlx + x~ , 77 = ^ln(x + /x^), ds 2 = dr 2 — r 2 dr] 2 . 
For any four-vector A p we have as A T = A T = (tA° — zA s )/t = (x + A~ + x~~A + )/t and 
A v = — r 2 ^ = zA° — tA 3 = x + A~ — x~ A + . This also applies to Dirac gamma matrices, 
giving 7 1 " = ^ e~ m0 ^ [i . We shall frequently separate Dirac spinors into eigenvectors of 7°7 3 , 
using the projection operators 

P± = 1(1 ± 7 V) = ^TV = -^7° = 5 7*7*, (1) 

satisfying P ± P ± = P ± , P ± P =F = 0, P + + P~ = 1. For momenta we use the transverse 
mass Up = pt 2 + m 2 and rapidity y = \ ln(p + /p~), giving the energy E p = u> p cosh y and the 
longitudinal momentum p z = oj p sinh y. 

2 General formulation of the problem in 3+1 dimensions; ini- 
tial conditions on the light cone 




Figure 1: Left: The background gauge field from the classical field model. In the Abelian case the 
field is a pure gauge also in region (3). Right: The two possible trajectories of (retarded) propagation 
of the fermion. The spinor starts from a negative energy state e lq ' x v(q) and gets a fc + -kick from the 
gluon field (1) on the a; + -axis and then a fc~-kick from the gluon field (2) on the cc _ -axis; or in the 
other order. Finally the spinor is projected onto positive energy states e~ lp ' x u(p) inside the future 
light cone. The label U^U^) refers only to the Abelian case, when the pure gauge field inside the 
future light cone is given by U(vJJi2) — f(2)^(i)- 
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We want to calculate the number of quark-antiquark pairs produced by the classical color 
fields in the model developed in Refs. 0IBj and solved numerically in Refs. In [T7| . 

it was shown that the average number of produced pairs can be expressed as follows: 

(nqq) = J ^f2E ( Qin 6 out(P)Wp) in ) 

/P ' /q u{p)T R {p,-q)v(q) , (2) 



(2tt) 3 2E p J (27r) 3 2£, 

where T R (p, —q) is the amputated retarded propagator of the quark in the external field, — q 
and p being respectively the incoming and outgoing momenta. Note that this quantity is 
distinct from the cross-section to produce one pair, which would require the time-ordered 
propagator of the quark. The difference between the two cases is explained in more detail in 
|17j . For computational purposes, it is in fact more convenient to write the previous formula 
in coordinate space. This can be achieved by using the standard machinery of reduction 
formulas, and one obtains: 

u(p)TJp,-q)v(q) = lim I d 3 ^e iiE " x °' p - x) u\p)^Jx° ,x), (3) 

x -^+oo J 

where ijjq^x ,^.) is the solution of the Dirac equation in the presence of the external field, 
with a negative energy free spinor as its initial boundary value: 

V> q (x° -» -oo, x) = e^^-^viq). (4) 

This formula can be trivially modified in order to use the proper time r instead of x°. 
Moreover, although it tells us that the on-shell pair production amplitude is obtained only in 
the limit of infinite time, one can also consider an extension of this formula where the limit 
is not taken, thereby defining a "time dependent pair-production amplitude" which could be 
used in order to probe the typical time-scale necessary to produce the fermions. We shall 
give several examples of this below (see Fig. Eland Fig. EJ). 

Note that in this formalism, one neglects the backreaction of the produced fermions on 
the color field. It is therefore a good approximation only if the fermions do not outnumber 
the gluons. In an Abelian theory the calculation only involves pure gauge fields and the 
pair production amplitude (see Eq. below) can be calculated analytically 4 , as shown in 
|16j . Furthermore, in an Abelian theory, the square of this generalized amplitude is, in fact, 
time-independent, suggesting that all the pairs are produced instantaneously at the collision 
time (for a collision at infinite energy). 

In the non- Abelian case the gauge fields are known analytically up to the future light cone 
(t = 0) and numerically for r > 0. Thus we will have to solve the Dirac equation numerically 

4 However, the evaluation of this amplitude is complicated by some infrared singularities that need to be 
properly regularized |24|. 
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for that region, with the initial condition at r = given by essentially the same calculation 
as in the Abelian case. 

Let us first review the classical field model of [7JIH1- Then we shall repeat the calculation of 
fermion-antifermion production from Abelian Weizsacker- Williams fields from Ref. |16j and 
use it to find the initial conditions at r = for our numerical calculation in the non- Abelian 
case. 

Let us assume we have two nuclei moving along the light cone, corresponding to a current 

= ^+£(x-)p (1) (x T ) +«(* + )p (2 )(x T ). (5) 

The two colour charge densities p( m )(xr) are, independently for the two nuclei, drawn from 
a random ensemble, which in the original McLerran-Venugopalan model is taken to be Gaus- 
sian: 

(pU^T)p\ m) (yT))=9 2 ^ m) S ab 5 2 ( XT -y T ), m = l,2, (6) 

where fi is a parameter describing the transverse density of color charges and can be related, 
up to a logarithmic uncertainty, to the saturation scale Q s |25| . More generally the charge 
distribution in Eq. © is not known, but its evolution when probing smaller Feynman x values 
in the nuclei can be calculated from the JIMWLK equation, see e.g. |26j . 

One first calculates in the light cone gauge {A + = for the nucleus moving in the +z 
direction, and A~ = for the nucleus moving in the — z direction) the pure gauge fields 
corresponding two the two nuclei: 

A U)M = V (m) (x T )a j C/ ( t m) (x T ), m = 1,2. (7) 

These depend on the Wilson lines U^(xt) given by 

U(rn)(x T ) = exp ^-ifi^4(x r )^ . (8) 

In a temporal gauge 5 A T = the initial condition at r = for the color fields Ay(r, xy) and 
A v (t, kt) is given by these pure gauge fields corresponding to the two nuclei: 

^(0,x T ) = 4 ) (x T )+^ (2) (x r ), 

A"(0,x r ) = |[4 ) (x r ),Aj 2) (x r )]. (9) 
One then solves the equations of motion 

[D fl ,F^} = (10) 

5 The A T — gauge coincides with the light-cone gauge A + = (resp. A~ — 0) if x + = (resp. x~ =0). 
This is why we can use the gauge field of the nuclei before the collision, in two different light-cone gauges, as 
an initial condition at r = for the gauge field in the A T = gauge. 
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using these initial conditions to find the fields at later times r > 0. In this gauge it is easy to 
find the Hamiltonian and thus the energy of a given field configuration. Additionally, fixing 
the Coulomb gauge in the transverse plane, • = 0, one can also define a multiplicity 
corresponding to the classical fields. 

Let us then turn to solving the Dirac equation in the Abelian case, following ^B]. In |16j 
it is solved separately in different gauges, covariant gauge (called singular gauge in ^H!) and 
light-cone gauge. In [2U] the covariant gauge is used. Here we use the light-cone gauge, 
because it is the same gauge that was used in solving the Yang-Mills equations up to the 
r = light cone. 

Following Eqs. © and (0} we start with a negative energy plane wave for x^ < : 



V>(4) (a 



iq-x 



v(q). 



Ill) 



The boundary conditions when crossing the light cones can be derived from the following 
argument. Since j + P~ = 7~ P + = 0, the Dirac equation only involves terms like d~P~ip 
and <9+-P + V) but not d-P + ip and d+P~ip. If there were a discontinuity in P~ip on the x~ = 
light cone, the derivative term would give a delta function, with no other term to compensate 
it, which is not possible. A discontinuity in P + ip, on the other hand, is possible because 
it can be compensated by the #(x~)-discontinuity in the gauge field. Thus the boundary 
condition is that on the x^ = light cone i/j^ is continuous. Using this boundary condition 
one can find the solutions in the regions (1) (x~ > > x + ) and (2) (x + > > x~): 



^(1)0*0 = ^(1)( x t) 



d 2 k- 

(27T) 



qrje 



x exp iq x 



2q- 



p- + p + i 



7r • k T 



m 



V2q- 



v(q) 



V>( 2 ) {x) = 17( 2 ) (x r ) 



d 2 k- 



(27T) 



q T Je 



-iky -Xy 



x exp 



2q 



-x + + iq + x 



p + + p~i 



o7r 



m 



v(q) 



(12) 



where oj\ 



k^ 2 + m 2 and we have defined the Fourier transforms as: 



d 2 y T e 



(13) 



Next we must continue to the forward light cone. The solution that matches to Eq. (|12l) 
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on the light cone is 



^(3)(s 



[/ (1) (x T )LT (2) (x T ) 



d 2 pr d 2 ky 
(2vr) 2 (2vr) 2 



car 



x exp 



P- 



2iri 



i | ^(2) (PT - kr)^ (k T - q T ) 

7t -kr 



P" 1 



2cr 



ie 



2^ 



+ p - 7 o7t-Pt 
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+ LT (1) (x T )[/(2)(x T 



y/2p+ - ie 
d 2 pr d 2 k T 



V2q- 



v{q) 



x cxp 



p- 




(14) 



Here the first term corresponds to a situation in which the positron state q first hits the 
nucleus 1 moving in the x + direction, propagates over region 1, meets the nucleus 2 and 
propagates into region 3 (the branch on the left in Fig. ^) . For QED the order of the Ur m \ 
matrices is irrelevant, not so for QCD. The p^-integrals in Eq. (|14jl can be performed to turn 
the expression into a sum of Bessel functions of the kind we shall encounter in Sec. H3 but we 
will not write down this complicated expression here. 

To find the matrix element for pair production, one has to project the spinor (|14|) to a posi- 
tive energy state e~ tp ' x u{p). Removing the product U( 2 )U(i), which is a gauge transformation 
to Coulomb gauge, one gets the Abelian theory result of |16j : 



M(p,q) = iV2 



d 2 k T [ UL(-PT ~ k T )uUk T - q T ) 



+ 



Lo q uj p eyp yi + uj\ 



T 



m) v(q) 



U^i-pr -k T )u} 2) (k T - cit) 



u\p)"l + (7 T • k T - m)v(q) >. (15) 



Kinematically, the terms in (|15j) correspond to the process k\ + k 2 — > p + q with k\ = 
{p + + q + ,0,'k.T + Pt), k,2 = (0,p~ + q~,q_r — kr) (see Fig. EJ). The two terms are the t- 
and li-channel propagator pole terms in the Feynman diagram corresponding to Fig. |5J with 



m~ 



m 



(ki —p) 2 = oJ q u! p e yq Vp + u>l and m 2 



m 



(ki - q) 2 = 0J q u) p e y p y " + u)\. 



We now wish to take the Abelian solution for r > 0, Eq. (|14j) . and write it in a form suitable 
for determining the initial condition as r — > 0. What is crucial here is the choice of the other 
variable, kept fixed when taking this limit. The choices are t], z, x . To obtain the correct 
result we must have a dimensionful longitudinal variable, such as z or to parametrise 
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ki 



P 



p 



+ 



q 



ki 



q 



Figure 2: Diagrams contributing to the lowest order pair production amplitude in QED. The incoming 
photons are quasi-real, with k\ = (p + + q + ,0, ky + Pt) and fc 2 = (0,p~~ + q~ , qr — kr). 



the r = surface. This is because one must be able to represent longitudinal momentum 
scales, for example uie y<1 /ujq, in coordinate space. For r > the corresponding longitudinal 
coordinate could be constructed as Te _r? , but for r = this is not possible. To enable a 
symmetric treatment of both branches in Fig. ^ we choose z as the longitudinal variable, 
with = (V T 2 + z 2 rfc z)/y/2 = (\z\ ± z)/\J~2 at r = 0. After a rather lengthy computation, 
the r — ► limit of the wave function Eq. (|14|) can be written as 



^(3)0" = 0,Z,X T ) 



d 2 k 



(27T) 



P + ^ r [/ (1) (x T )f/ ( t 1) (k T )exp 



2lj„ 



t — 1 



+P- 7 U [t 7T ■ D T - m] U {1) (xt)U} n(k T )-2 



cxp 



-.'A; 



+P" — C7 (2) (x r )f/f 2) (k r ) exp 



^j +g e-^(|z|+z) 



t . , 1 



+P+ 7 U [t 7 T • Dt - m] C/ (2) (x T )C/ ( T (k T ) — 



k+g 



cxp 



^fc+ q e^(|g| -z) 
2cj„ 



.^ +(? e-^(kl+^) 



2uj n 



X7° (7t • ( k T + Or) - "*) «(?), 
= Vt + i5Ar(0,xr) and uj\ 



(16) 



where D T = Vt + igA r (0,x T ) and = (k T + qr) + rrr. Note that whereas Eq. (fTljl 
involved products of £/m and J7(2), implicitly assuming that they commute, the terms in 
Eq. l(Trj|) only contain either Um or t/> 2 ) and thus Eq. (|Td)) can be directly generalised to the 
non-Abelian theory. The products of Um and Up) show up in the gauge field in the covariant 
derivative Dy, which depends on both Un\ and Um\ by Eqs. (JTJ),©. 

Now we have the initial conditions for a numerical solution of the Dirac equation in 3+1 di- 
mensions. To proceed further, one will have to generate the random SU(3) matrices £7( m ) (x^), 
compute the color fields A v (t, xy), Aj(r, xy), compute the spinor i/)rg\(T, z, xy) from the Dirac 
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equation using these color fields and the initial condition Eq. (jlfij) . project to a positive energy 
state e~ tp ' x u{p) and, finally, integrate the square of the amplitude so obtained over momenta. 
As this is a rather involved operation, we shall first simplify the problem by neglecting the 
transverse dimension. Computation of qq production in this 1+1 dimensional toy model 
permits us to test numerical aspects of the solution and the projection to final states. We 
shall return to the 3+1 dimensional case in future work. 

As a first step, setting C/( m )(xr) = 1, Ur m \{k.T) = (2vr) 2 5 2 (kj-) , in Eq. (fTB|) gives 



V>(3)( r = 0,z,x T )|{7 = i = e" 



Jq+x p+ _ ]J?_ e -yq / e iq + x _ -|\ ^0 p+ 

UJg 

+e iq-x+p- _ ™ e y q ( e iq-x+ -iWOp- v ( q ) (17) 

This also illustrates the structure of Eq. (|16|) : the left branch (first line) has for z < (on 
the x~-axis) firstly a component P + v(q) moving in the +z direction, but one also needs a 
component ~ j°P + v (q) to satisfy the Dirac equation. This is worked out explicitly below in 
Eq. (|29|) , For z > (on the x + axis) only e _iqT ' XT P + f (g) moving in +z direction remains. 
The right branch behaves symmetrically. 

3 Free Dirac equation in 1+1 dimensions; analytic 

Let us first define our spinor conventions and study the solutions of the free Dirac equation 
in 1+1 dimensions. 

In 1+1 dimensions, given fermions of mass m, we can parametrise an on-shell momen- 
tum vector by just the rapidity y: (E,p z ) = m(coshy, sinhy). A free wave is then e ip ' x = 

^irriT cosh(y— r]) 

Let us choose for the Dirac matrices a representation where 7°7 3 is diagonal: 

^-(r;)- "-.«-(;_!)■ m 

In this basis the projection operators defined in Eq. Q are simply: 

Hoc)' 

and we denote the two components by 

*=(£)■ <20) 

The Dirac equation: 

{vfdp -m)ip = (21) 
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has plane wave solutions corresponding to positive and negative energy: 



^ {+) {x) = e-^ x u{y) 4> H (x) = ^v(y)- (22) 

Using the explicit forms of the Dirac matrices we can easily see that 



u{y) = Vm[_i y j, v(y) = y/mi -i y J- (23) 

The solutions have a Lorentz-invariant normalisation with u{y)u{y) = v{y)v{y) = 2m and 
u(y)v(y) = v(y)u{y) = 0. The quantity we will be interested in, however, is the particle 
number density on a constant proper time surface, which is the r-component of a Lorentz 
vector. For both the positive and negative energy solutions it is given by u(y)'j T u(y) = 
v(y) , y T v(y) = 2mcosh(y — rj). The cross term u(y)^ T v(y) = 2msinh(y — rj) is not zero, but 
vanishes by symmetry when integrated over the longitudinal coordinate rj. 

Writing the free Dirac equation in terms of ift = P ip, the eigenvectors of 7°7 3 , we have 

a T + ^ ± = me v, (24) 

or, squaring to get a second order equation, 

« 2 1 

^ = 0. (25) 



d T i + U T -%+m 2 

T T z 



With an ansatz ip^ 1 (r, 77) = e nr, ip^ (r) this reduces to the Bessel equation. The solutions that 
are separable and finite for r = are of the form e nri J n (mT). 

To find the right linear combination of the separable solutions we have to look at the initial 
condition. This can be done by looking at the full initial condition, Eq. (|16|). and removing 
all the transverse degrees of freedom. One first takes U^(xt) = 1, which leads to Eq. (fTTj) . 
and further sets = and thus reduces uj q to m. Using the 2d representation of 7 , 7 3 and 
v(q) given above, Eq. (fTTj) then becomes 6 

/ e y/2 e imey(\z\-z)/2 \ f e y/2 Lime-v (\z\+z)/2 _ A \ 

*<J = 0,Z)=V^[^ _ e -y/2^ e ime"Qz\-z)/2 _ 1} j + ^-y^m^z^ ' ) • 

(26) 

where the first(second) term is the left(right) branch in Fig. ^ From this we find that the 
initial condition for, for example, the upper component of the left hand branch, must behave, 
up to a sign, for r — > as: 

^+( r 0,77) = -^ e y/ 2 e^ mTeV ~" = - y / m - e y / 2 e i > eV ( V:fT +* 2 - z \ (27) 

6 We will actually change the sign of this expression; this would correspond to changing the sign of v(q) in 
Eq. i'2'M . which is just a convention. 
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Such a solution can be constructed as a sum of Bessel function modes as 



mr); 



(28) 



n=0 



using J n (mr) — > (mr/2) n /n! for r — > and summing over n gives Eq. (|27|) . Using Eq. (|24j) 
we see that the other component of the spinor is: 



mr) 



(29) 



n=l 



which is exactly the same as the lower component of the first term in Eq. (|26|). showing 
the consistency of the approach. We show in Appendix [B] that when this wave function is 
projected on positive energy states, one obtains the same result (±i/ cosh ^ (y — y') for the 
two branches) for the amplitude as from evaluating the Feynman diagrams in Fig. |21 or by 
specialising the Abelian amplitude in Eq. (fT3)) to 2d. 



4 Free Dirac equation in 1+1 dimensions; numerical 

Let us now formulate the discretised solution of the Dirac equation using the coordinates 
r, z. The advantage of this coordinate system is that it allows a simultaneous and symmetric 
treatment of both two branches. On the other hand, at r = 0, yplx^ 1 = \z\ ± z is not a 
continuous function of z and there are corresponding discontinuities in tp. The free Dirac 
equation in this coordinate system is 



d T ip 



± 



\Jt 2 + z 2 ± z 



-3//2 



exp i 



The initial condition for the left hand branch is: 

V,+ (r = 0,z) : 

T/T(T = 0,Z) : 

and for the other one: 

1P + (T = 0,Z) : 

iI>-(t = 0,z) -- 



i 2 I me y 
e y/ exp ( i [\z\ 



me" 



. .me y . 
1 — exp ( i — - — (\z\ + z) 



-3//2 



exp i 



.me 



\z\+z) 



(30) 

(31) 
(32) 

(33) 
(34) 
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Because the z-dependent coefficients in the equation would make any explicit scheme un- 
stable, we discretise Eq. ()3()[) implicitly 7 : 

1 



2dr 



[ip* (r + dr, z) - V ± (r - dr, , 



Vt 2 + z 2 ± z 



ip (r + dr, z + dz) + ip (r — dr, z + dz) 



Ardz 

-ip^ (r + dr, z — dz) — ip (r — dr, z — dz) 



■ Vt 2 + z 2 ±z 

vm W \ T i z ) ■ \yv) 



Now the ±-components are stored at different timesteps: ip + (r); ip~(r + dr). This saves 
memory compared to storing the spinor at two timesteps, while the discretisation is still 
second order accurate in dr. It seems that it is critical for the stability of the algorithm to 
also discretise the endpoints to second order accuracy in dz. Thus, when in Eq. 1)35(1 we have 
used the centered difference 

f'(z)^^-[f(z + dz)-f(z-dz)], (36) 

for the points inside the lattice, for the edges of the lattice Eq. must be modified to use 
a one-sided second order accurate difference: 



2f(z + dz)-±f(z + 2dz)-p(z) 



(37) 



This prescription could be described as a free boundary condition. Note that it would be 
quite unphysical to impose periodic boundary conditions in the z-direction. Technically this 
shows up e.g. in the fact that the coefficient in front of the <9 z -term in Eq. (|30|) would be 
discontinuous at such a periodic boundary. 

Eq. (|35|) forms a system of linear equations for ip^ (r + dr) in terms of the known values 
-0 =t (r — dr) and ip^(r). The system is almost tridiagonal and can be efficiently solved using 
LU-decomposition, leading to an algorithm which is slower than the corresponding explicit 
discretisation only by a constant factor. 

After having solved numerically the Dirac equation to find the spinor at some finite proper 
time r, we must project out the positive energy part of the wavefunction to find the amplitude 
for production of a pair at rapidities y,y' . This is given by 

dz 



M(y',y) 
where 



Vt 2 + z 2 



u(y')exjp(im \J t 2 + z 2 cosh y — z sinh y ]7 r ^(r,z 



(38) 



7' = 7°e vi ! 3 = jO (gQgh^ _ j 3 sinh 7]) = "f" | 7"7 

V r V 

and r/V r 2 + z 2 is the Jacobian. 

7 The discretisation of a partial differential equation is said to be "implicit" when the time derivative of the 
unknown function at a given timestep depends on the unknown function at the next timestep. 



\/r 2 + z s 



0„3 ' 



(39) 
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Figure 3: Numerically calculated real and imaginary parts of the free amplitude M for one branch 
shown for different lattice sizes in the z-direction. Also shown is the analytical value 1/ cosh(y — y 1 ) 
of the imaginary part. The analytical value of the real part is zero. 

In practice the integral, or sum in the discrete case, is an oscillatory function for large z. 
In an analytical integration these oscillations average to zero, but in a numerical calculation 
with a finite extent in the z-direction this is harder to achieve. We have used two different 
techniques to treat this problem. One is to calculate the integral (|3"%|) for different upper 
and lower limits ±z max , and then take an average of the values thus obtained over a range in 
z max that contains several periods of oscillation. The other technique is to use wave packets 
which are localised in the ^-direction to slightly less than the extent of the system in the in- 
direction. This introduces an uncertainty to the momentum in the z-direction or equivalently 
the rapidity y, but the uncertainty is of the same order of magnitude as the infrared cutoff 
from the size of the z-lattice. 

The range of rapidity y that can be reached in this coordinate system is limited by two 
things. The finite lattice spacing in the z-direction gives an ultraviolet cutoff for p z , implying 
that we must have sinhy < l/(mdz). Our method also requires that the values of rj of 
the order of the momentum space rapidities studied be covered by the finite lattice in the 
z-direction. This translates into the requirement sinhy < z max /r. Fig. Olgives an idea of how 
close our results are to the analytically known result (for one branch) M(y' ,y) = ij cosh \{y — 
y') (see Appendix |B|). 
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5 Dirac equation in 1+1 dimensions; external field 



First let us note that this is not a calculation of pair production from the Weizsacker- Williams 
fields of two currents on the light cone, because such a thing does not exist in two spacetime 
dimensions. Assuming an external current 

J + = e\5(x~), J~ = e2<5(x + ), 

one can namely solve d^F^ = J u for the only component E = to be E = — e\Q(x~) + 
e2©(x + ). There is just a constant electric field off the light cone. Instead, our purpose is to 
impose by hand an external field to test our numerical method and to model the projection 
of the real 3+1 dimensional physics on the longitudinal dimension. The parameters of the 
calculation will be m, the mass parameter of the Dirac equation, Q s , another mass parameter 
describing the proper time variation of the external field and c, a dimensionless parameter 
describing the strength of the external field. To effectively describe the omitted transverse 
momentum effects one might take m to be not very different from Q s . For c -C 1 one is in 
the weak field domain and analytic results can be obtained. 

With the gauge choice A T = and an external gauge field A ti (t) the Dirac equation (|3l7|) 
becomes: 

dr1 p± = Vr2 + z2±z ( Td ^± _ im ^\ T iiMll^. (4 ) 

T T 

We shall study this for various choices of A v (t). The first choice, motivated by the pertur- 
bative solution of the Yang- Mills equations [Jj is 

gA r ,(r) = cQ s TJ 1 (Q s T). (41) 

This form is special in that the Fourier transformed fields corresponding to (J41|) can be given 
analytically: 

^ ±(t+ , fc - ) = ± ^^___i__. (42) 

Note that the pole structure is dictated by the requirement that A^ ~ 9{x Jr )9{x^). Using 
this explicit form for the field in (|41|) and the spinors in (|23|) one can write the lowest order 
perturbative result (corresponding to diagram (a) in Fig. EJ) for the amplitude: 

M = -igu{p)4{p + q)v(q) (43) 

in the form 

M (A = - ') = ^ 2cQ s 2 

{ V V V ' cosh(i Ay) [2m 2 (l + cosh Ay) - Q 2 S + ie] cosh(i Ay) [s - Q 2 S + ie] ' 

(44) 

Numerical results for the field 1)41(1 are shown in Fig. [1] at a fixed large time r = N T dr, the 
time dependence is studied in Fig. El In Fig. 0] the left panel shows the amplitude for weak 
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Figure 4: Absolute value of the quark pair production amplitude for different values of the oscillation 
scale Q s . Left: weak fields (cQ s — 0.05m), the peaks are at the location given by Eq. Q45J1. Right: 
strong fields (cQ s — m) with the same values of Q s . Peaks near the "threshold" Q s — 2m are shifted. 

fields and the right one for strong fields. For weak fields, c<l, one expects the numerical 
result to coincide with first order perturbation theory, as given by (j44j) . which has a peak at 
s = Q 2 or at 

cosh^ = ^. (45) 

For Q s < 2m this equation has no solution for Ay and, in fact, the numerical result is very 
small. Precisely at "threshold", Q s = 2m, there is a very strong single peak at Ay = 0, 
the quark and antiquark emerge at rest relative to each other. For Q s > 2m there are two 
peaks corresponding to the two signs of solutions of These are well reproduced by 

the numerical calculation. Due to the finite time the peak is not a delta function, but is 
broadened. Physically, the amplitude peaks at pair invariant mass = Q s and, in 1+ld, the 
only way to give the pair this invariant mass is to separate them in rapidity. In 3+ld the 
situation is quite different, since then 

s= (p + q) 2 = 2uJ p LJq cosh Ay + 2m 2 - 2p T ■ qr (46) 

and the pair invariant mass can even be dominated by transverse momenta. 

For stronger fields the numerical calculation sums over all orders in the external field (all 
diagrams in Fig. and one does not necessarily expect any peak structure. However, peaks 
still appear (right panel of Fig. 0]), although the location of the peaks is shifted, especially 
near Q s = 2m. 

Fig. shows the dependence on the physical time r = N T dr, at which the projection to 
the final state is done. One sees that the height of the peaks increases essentially linearly 
in time, M pea k ~ N T , while their widths shrink somewhat. The area of each of the peaks, 
J dy|M|, is approximately independent of N T , showing that the resulting rapidity distribution 
at asymptotic times contains two delta function peaks in Ay. The right panel shows that the 
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Figure 5: Diagrams contributing to the amplitude in the 1+1-dimensional toy model 
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Figure 6: Left: Absolute value of the amplitude for Q s = 2.05m and cQ s = 0.5m when the projection 
is made at different physical times N T dr at a fixed dr = 0.05/m. Right: The number of produced pairs 
per unit rapidity, J dAy\M | 2 , as a function of time. Note that Q s > 2m, so that the delta-function 
peak in Eq. (|44H dominates. 



integral J dy\M\ 2 , giving the number of produced pairs, grows approximately linearly with 
r. This monotonic increase is due to the fact that the ansatz of Eq. (j41j) behaves like ~ yfr 
at large r. 

As a second example we shall consider a non-oscillatory and exponentially decaying field 



gA v (T) = cQ s re QbT . 



(47) 



This is actually simply reproduced from the first ansatz of Eq. (|41|) by taking a superposition 
of Bessel functions ujtJi(ujt) with different frequencies u, which "washes out" the peaks at 
s = to 2 in (|44j) . The appropriate weight factor is given by the relation 



9M T ) 



r i d0O- 



Q s uj 



r wTJi(wr) = cQ s re QsT . 



(48) 



Integrating over the matrix element (|44j) (with Q s — > uj) one finds the perturbative matrix 
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element 
M(Ay) 



cQ s s 



cosh(±Ay)(s + Q s 2 ) 3 / 2 



m + 2 



(49) 

Now there is no peak at s = Q 2 and, when plotted as a function of Ay for various Q s /m, 
M(Ay) is a Gaussian-like curve centered around Ay = 0. 
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Figure 7: Left: Numerically computed amplitude for the strong (c = ±1) exponentially decaying field 
in Eq. (|48|l . The curves labeled "theory" are the weak field analytical result from Eq. I|49|l . Changing 
c — * —c reflects the curves around Ay = 0. Right: The number of produced pairs per unit rapidity, 
J dAy\M\ 2 , as a function of time for the exponentially decaying field. 



Numerical results for the amplitude for Q s ~ m and a strong field c = ±1 are shown in the 
left panel of Fig. [7| The general form of the rapidity dependence and the normalisation agree 
quite well with the weak field formula (|49|) , but the numerical curves are not centered exactly 
around Ay = 0. This asymmetry under Ay — > —Ay is simply due to the fact that under 
parity A v — ► —A v and thus the ansatz (|47|) (as well as any non-zero A„) breaks parity. To 
check that this is a real physical effect the numerical calculations have been performed with 
both c = 1 and c = — 1. In the Dirac equation this corresponds to z — > — z and Ay — > — Ay 
and the curves obtained agree with this. The right panel of Fig. shows the number of pairs 
produced, J dAy|M| 2 . Initially it rises ~ Q s t, then after a few oscillations, caused by the 
coupling of the two frequency scales Q s and m, reaches a constant multiplicity level. Q s being 
a timescale of damping, the oscillation frequency is given by m. The larger Q B /m, the larger 
is the multiplicity. This constancy is due to the fact that the external field vanishes in a time 

~ i/Qs- 
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6 Numerical tests 



Because of the delicate nature of the time-longitudinal dynamics, we have performed various 
tests of the numerics. In Fig. El we study the effect of taking different (physical) sizes for the 
lattice in the z-direction, or z max , using the Bessel function external field l|41|). Especially 
when the projection is done at a larger time the amplitude depends somewhat on the size 
of the lattice. The difference in the multiplicity J dy|M| 2 is, however, quite small. For the 
exponentially decaying field l)47|) the lattice size effect is much smaller, almost unobservable 
on a plot like Fig. |H1 The effect of a finite z max on the integral J dy\M\ 2 is larger for the 
exponentially decaying field, because of the contribution from larger values of |Ay| that are 
unaccessible for a small z max (see the discussion at the end of Sec. HJ. 




) , , , , 1 , , , , 1 , , , , 1 , , , , 1 II 1 1 1 1 T 1 1 1 1 1 ' 1 1 
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Figure 8: Absolute value of the amplitude for different sizes of the lattice size in the z-direction. The 
external field is a Bessel function Eq. (|41|) with a frequency slightly above the resonance condition, 
Q s = 2.05m. Left: c:Q s = 0.615m and N T dr = 60/ra. Right: cQ s = 0.5m and N T dr = 100/m. There 
is some dependence on the lattice size, and the dependence is larger when the projection is done at a 
later time (right panel). 

The left panel of Fig. |2l shows the dependence on the timestep used. In the right panel of 
Fig. E3 we choose different rapidities y for the antiquark and compute the distribution in the 
quark rapidity y' . The outcome is always a function of y — y', which shows that our numerical 
method preserves the boost invariance of the result to a good accuracy. In Fig.EHwe explore 
the accuracy that can be reached with lattices small enough to make a full 3+1-dimensional 
computation realistic. Although the computational requirements of a 3+ ld-simulation are 
quite hard, we believe that with a careful choice of discretisation parameters it is possible to 
extract some physical results from a full 3+ Id numerical calculation. 
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Figure 9: Left: Absolute value of the amplitude using different timesteps at fixed physical size 
N T dr = 400/m. Right: amplitude for different values of the antiquark rapidity y; to check that the 
numerical calculation reproduces the boost invariance of the solution. Both plots have an oscillating 
field with Q s = 2.05m and cQ s = 0.5m but different lattice sizes in the ^-direction. 
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Figure 10: Absolute value of the amplitude for smaller lattice sizes in the z-direction. Left: Q s below 
resonance (Q s = m and cQ s = m). Right: Q s above resonance (Q s = 2.05m and cQ s = 0.5m). Values 
of z in units of [1/m]. 

7 Conclusions 

In this paper, we have set up the framework for computing qq pair production in ultrarela- 
tivistic heavy ion collisions in the classical field model with an ensemble of quantum initial 
conditions. This is an important theoretical problem, especially in view of nonperturbative 
chemical thermalisation, for which one needs a sufficient number of qq pairs from the domi- 
nantly gluonic initial state. However, the calculation is technically complicated, involving the 
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numerical solution of both the gauge field equations of motion and the Dirac equation in the 
background gauge field. Even the formulation of the initial condition proves to be nontrivial 
since the natural variables r, rj cannot in the limit r — > give a dimensionful longitudinal 
variable, which one needs for longitudinal Fourier transforms. In view of this, we have in 
this paper limited ourselves to giving only the initial condition for the full 3+ld problem but 
considered in numerical detail only a 1+ld version of the model obtained after truncation of 
the transverse dynamics. In this model both rapidity distributions and the total number of 
produced pairs were computed for two forms of the gluonic external field. 

The work carried out here has solved the conceptually most complicated part of the full 
3+ld problem, the formulation of the initial condition and the treatment of the longitudinal 
dimension together with the proper time. The inclusion of transverse dynamics is compu- 
tationally demanding but otherwise straightforward. When that part is completed, one will 
be in a position to make meaningful statements about the nonperturbative production of 
qq pairs in heavy ion collisions. 

Future tasks include a full 3+1-dimensional treatment of also the gauge field equations of 
motion and ultimately also including feedback from the qq sector to gluons, i.e., formulating 
and solving coupled equations of motion. This will be a very challenging task. 
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A Dirac equation in curved coordinates 

Let us denote the flat coordinates t,z or x°,x 3 by Latin indices a, b,... and the curved 
ones r,rj by Greek ones: fj,,u,.... The flat metric is rj^ = diag(l, — 1) and the curved 
one g^v = diag(l,— t 2 ), g^ v = diag(l, — 1/r 2 ). The nonzero Christoffel symbols for the r, r\ 
coordinates are [2Z] VL„ = r and = T^ T = 1/r. 

Given some representation for the usual 7-matrices in flat space, j a , one can express the 
7-matrices in curved coordinates as 7^ = ea"f a . The zweibein ea relates the flat metric to the 
curved one by g^ u = e^e^r^ and, conversely, r\ a ^ = e^e^g^. There is no unique choice for 
the e^, reflecting the fact that there are different ways one can attach a flat tangent space to 
each point in spacetime. We have mostly used the natural intuitive choice for the zweibein, 
namely e° = d^x a . But in the r, ^-coordinate system there is also another natural choice, 
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namely to take e° = 1, ef. = r, so that 7^ do not depend on the coordinates. To preserve the 
local Lorentz invariance of the Dirac equation, one must introduce a spin connectAon[W\: 

r, = 1^,^(8^ + ^). (50) 

In this case the spin connection has only one nonvanishing component: 

IN, = ^7°7 3 - (51) 
The free Dirac equation [i^id^ + r„) — m] i/j = in this case becomes: 



m 



ip = 0. (52) 



Here we have introduced the spinor i/j defined with this choice of the zweibein. It is related 

10 3 

to the usual flat space spinor by ip = e - ^ 7 7 ip. The plane wave solutions (122H and (|23j) now 
have a form that makes boost invariance manifest: 



*i±)W = V^e^co sHy - v) ( ± e ^ ) ) • (53) 



B Amplitude 



We found that the solution of the free Dirac equation for the left branch in Fig. ^ in the 
future light cone is (Eqs. (J2HJ) and (|29|0 



00 



= -vW/^^-TJ^mr) (54) 

n=0 

00 

</r = ^R e - y ^^2{ie y ~^) n J n {mT). (55) 

n=l 

To project to a positive energy state with rapidity y', we calculate the amplitude 

M 1 (y',y) =r J dr/^V^^VV'M). (56) 
Inserting the form (|54j) and using the standard integral 

dr} e K cosh^-j,)-^ = e - uy me iu *l 2 H^{Q (57) 
and the Wronskian relation 

(mr) J n+1 (mr) - fl$l (mr ) J n (mr ) = — (58) 
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one gets 



M 1 (y',y)= * . (59) 
cosh M -^ L - 



The other right branch has the initial condition 

oo 

= ^-^^(ie^^V^mr) (60) 

n=0 

oo 

^+ = -^/ 2 ^(ie"-») n J n ( m r) (61) 
n=l 

and leads to a contribution which exactly cancels that from the left branch: 

M 2 (y',y) = — l — T . (62) 
cosh 

These two contributions and the way they cancel are exactly the same that come from eval- 
uating, in 1+1 dimensions, the Feynman diagrams in Fig. 
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